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We calculate the ground-state properties of unpolarized two-component Fermi gas by the diffusion 
quantum Monte Carlo (DMC) methods. Using an extrapolation to the zero effective range of the 
attractive two-particle interaction, we find E/E^^^ to be 0.212(2), 0.407(2), 0.409(3) and 0.398(3) 
for 4, 14, 38 and 66 atoms, respectively. Our results indicate that the dependence of the total 
energy on the effective range is sizable and the extrapolation is therefore quite important. In order 
to test the quality of nodal surfaces and to estimate the impact of the fixed-node approximation we 
perform released-node DMC calculations for 4 and 14 atoms. Analysis of the released-node and the 
fixed-node results suggests that the main sources of the fixed-node errors are long-range correlations 
which are difficult to sample in the released-node approaches due to the fast growth of the bosonic 
noise. Besides energies, we evaluate the two-body density matrix and the condensate fraction. We 
find that the condensate fraction for the 66 atom system converges to 0.56(1) after the extrapolation 
to the zero interaction range. 



I. INTRODUCTION 

In recent years, the homogeneous Fermi gas with at- 
tractive interactions has been studied extensively both 
theoretically and experimentally due to the success in 
cooling atoms into ultracold dilute condensates [TH3]. 
By tuning the interaction strength through the Feshbach 
resonance, 4 - 7 the system can cross from the Bardeen- 
Cooper-Schrieffer (BCS) superfluid phase, where the 
s-wave scattering length is negative, to the Bose- 
Einstein condensate (BEC), where Os is positive. Since 
there is no symmetry change of the quantum state in- 
volved, the system exhibits the well-known BCS-BEC 
crossover. 

In the special case corresponding to the diverging scat- 
tering length, flg — >■ oo, the system is in a strongly inter- 
acting regime called the unitary limit. In this regime the 
intcrparticle spacing is the only relevant scale, and 
the rest of the quantities are universal and system inde- 
pendent. The total energy of this system can be conve- 
niently written as i? = ^iJfico, where i?frco is the energy 
of the non-interacting atomic gas and ^ is a system in- 
dependent parameter. Experimental measurements of ^ 
have been performed using ^Li and ^''K atoms by inves- 
tigating the expansion rate of the atomic cloud and the 
sound propagation in it[8l412]. Simultaneously, a number 
of theoretical and numerical estimations of ^ have been 
reported, including diffusion Monte Carlo (DMC) [TMT5] 
as well as path integral Monte Carlo, lattice simulations 
and analytical methods |18H27j . The resulting estimates 
fall between « 0.25-0.45 showing that the actual value 
has not been settled yet and is still of significant interest 
due to the universal nature of the unitary limit. 

One of the most interesting properties of the unitary 
gas is the robust presence of the pairing condensate which 
involves a large fraction of the system. The study of pair- 
ing effects is thus much more straightforward than, say, 
in superconducting materials, where only a sliver of the 



fermions around the Fermi level forms the condensate 
and the attractive interaction is much more complicated. 
The quantum Monte Carlo (QMC) methods have the ad- 
vantage that the condensate can be detected directly, by 
evaluating the off'-diagonal two-particle density matrix 
and by monitoring its behavior at large distances [TH1I18) . 

The goal of our study is twofold. First, any actual sim- 
ulation involves only a finite system, and the quantities 
relevant for the thermodynamic limit have to be obtained 
from appropriate extrapolations. The unitary system is 
not trivial in this respect, since pairing with infinite scat- 
tering length is described by a function with a slow fall-off 
at large distances. It is necessary to analyze the finite-size 
scaling of the quantities of interest and to test whether 
the actual limit of infinite dilution, or, equivalently, of 
point-like character of the interaction, has indeed been 
reached. Second, the impact of the fixed-node approxi- 
mation in the quantum Monte Carlo method is not very 
well understood for this system since there is nothing to 
compare with: so far the fixed-node formulation of the 
QMC methods appears to be the only approach that is 
able to provide an upper bound for the total energy. This 
has motivated us to probe the accuracy of the nodes by 
released-node QMC simulations and by improvements in 
the variational fiexibility of the employed wave functions. 

We have carried out calculations of the ground-state 
properties of the dilute unitary Fermi gas by the fixed- 
node DMC (FN-DMC)[5Hj method for 4, 14, 38 and 66 
atoms. By using a more versatile construction of the pair 
orbital in the BCS wave function and by extrapolating 
the effective range of the two-particle interactions i?eff 
to zero, we were able to obtain lower ^ than, for exam- 
ple, the previously reported DMC calculations in Ref. 
[TB] based on BCS trial functions. These results suggest 
that the extrapolation of i?cff is important, especially for 
smaller systems. In order to test the quality of the nodal 
surface of the BCS wave function, we have performed 
released-node DMC (RN-DMC)[ini EO] calculations for 
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4 and 14 atoms. This procedure has been carried out 
starting from two types of nodal constraints: from the 
BCS nodes and from the Hartree-Fock (HF) nodes. Our 
RN-DMC results indicate that the nodal corrections are 
driven mainly by long-range correlations which are diffi- 
cult to sample in the released-nodc framework due to the 
rapid growth of the bosonic noise. We have calculated 
also the two-body density matrix and the condensate 
fraction for the 66 atom system, and we have estimated 
the corrections from the effective-range extrapolation on 
these quantities. 



II. METHOD 

A. Hamiltonian 

We consider a two-component Fermi gas with Hamil- 
tonian 

N/2 N/2 

^ = -^Ev?-^Ev?'+E^(--), (1) 

2—1 i' — 1 

where N is the total number of atoms, i and i' correspond 
to the spin-up and spin-down atoms, and rai denotes the 
distance jr^ — ri/|. The atoms are located in a cubic box 
with the side L and we impose the periodic boundary 
conditions. The two-particle potential V^ru') is taken in 
the Poschl- Teller form 



V{ru') 



cosh {firw) 



(2) 



whose effective range is i?cfF — 2//i. The s-wave scatter- 
ing length as is infinite for all values of /i 7^ 0. 



B. Trial wave functions 

In the majority of our calculations we employ trial 
wave functions of the BCS form multiplied with the Jas- 
trow factor (BCS-Jastrow) as given by 



(3) 



where 



r N/2 



A 



det[0(r„r,O]. (4) 



Here A represents the antisymmetrization operator and 
(j){ri,rir) is the pair orbital. The vector R encom- 
passes all atomic coordinates and r^/ . Additionally, we 
have carried out a subset of calculations also with the 
Hartree-Fock-Jastrow (HF-Jastrow) trial functions, in 
which ^E" BCS is replaced with a product of two Slater de- 
terminants of one-particle orbitals (simple plane waves). 
The HF-Jastrow wave function reads as 



The pair orbital (j){ri, r^') in ^bcs is written as a linear 
combination of Gaussian functions 



l,rn,n— — l k 



where dk are expansion coefficients, and = (xi,yi,Zi) 
and Ti' = {xi> , , z^/ ) are coordinates of i and i' atoms 
inside the simulation box of linear size L. We choose suffi- 
ciently large exponents ak so that only the first neighbor 
shell of periodic images contributes to the sum, that is, 
the Gaussian functions are negligible at distances larger 
than 3L/2. The pair orbital is smooth with zero deriva- 
tive at the boundary of the simulation cell. The Jastrow 
factor J(R) is constructed in a similar way as the pair 
orbital 0(ri,ri') for both different spin atoms and same 
spin atoms. 

A typical trial wave function includes around 30 to 
40 variational parameters that are optimized by min- 
imizing a linear combination of the total energy and 
its variance[ni]. Although the Jastrow factor does not 
change the nodal surface, accurate description of the pair 
correlations makes the variational optimization much 
more efficient and robust. When the effective range of the 
potential approaches zero, more Gaussian functions with 
larger exponents ak are included in the Jastrow factor in 
order to keep the accuracy of the trial function consis- 
tently high. On the other hand, and somewhat surpris- 
ingly, we find that similar adjustment of the pair orbital 
with the changing effective range is relatively minor. 



C. Fixed-node and released-node DMC methods 

The DMC method projects out the ground state from 
a given trial function by means of an auxiliary evo- 
lution in the imaginary time, $(t) ~ exp(— r-ff)\E'T. By 
introducing importance sampling 28J with the aid of a 
guiding function ^Po, we can write an integral equation 
for $(R, r) in the form 

*g(R)$(R,t + At) = 

/ ^R'|f^G(R,R',AT)vI/G(R')'i>(R',r). (7) 

For small At, the propagator G(R, R', At) can be ap- 
proximated using the Trotter-Suzuki formula as 

*^^^G(R, R', At) « Go(R, R' + Atv(R'), At) 



X e 



-Ar[EL(R)+EL{R')-2ET]/2 



(8) 



(5) 



where v(R') = Vln|*G(R')l and Go(R,R',At) is the 
Green's function for non-interacting atoms that takes the 
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FIG. 1. The fixed-node energy for unpolarized unitary Fermi 
gas as a function of tlie interaction range Raa/vs witli linear 
extrapolation to R^g/rs — 0. The system sizes are 4, 14, 38 
and 66 atoms from the top left to the bottom right. The 
statistical error bars are smaller than the symbol size. 



form of the diffusion kernel. The so-called local energy 
El is given by 



*g(R.) ■ 



(9) 



The product ^'c'l) is represented by a set of sam- 
ples (also referred to as walkers) and this ensemble is 
evolved with the aid of a stochastic process simulat- 
ing Eqs. Q and ([s]). In the fixed-node method we set 
^g(R-) = ^t(R) and the fixed- node condition is im- 
posed by enforcing the sampling points to obey 



*g(R.)*(R,t-) > 



(10) 



at all times. In the limit of long r the solution con- 
verges towards the lowest-energy state consistent with 
the boundary conditions given by the fixed nodes. 

In the RN-DMC method the guiding function has 
bosonic symmetry and its square should be close to the 
square of the fermionic ground state. We have used guid- 
ing functions in the form ^301 15^ 



vI'g(R) = V*t(R) + "(^i 



(11) 



where {'^t) is the average value of ^'|.(Ro) over all con- 
figurations, and Rq are the walker positions right after 
the nodal release. The tunable parameter a controls the 
rate of walkers passing through the nodal region. The 
guiding function is non-negative everywhere and there- 
fore the stochastic process propagates a mix of bosonic 



and fermionic states. The fermionic component is fil- 
tered out by reweighting with the factor so that 
the fermionic-state energy is given by 

/dR<i>o(R)*G(R)|4|i^PW 
($o|ff|*T) = ^ ^^^^ Z". . (12) 

/dR$„(R)*G(R)t|S 

where $0(1^) denotes the exact fermionic ground state. 

Since this method is exponentially demanding both in 
the projection time and in the number of atoms, it is im- 
portant to choose a so that the statistical information is 
recovered as quickly as possible. If a is too large the fluc- 
tuations from the poor importance sampling overwhelm 
any useful signal very rapidly. On the other hand, a 
too small value can bias the results. Since it is difficult 
to reach reliable error bars in this type of calculations, 
we have used the method mostly to identify the onset 
and the amplitude of the energy decrease during the pro- 
jection period when the stochastic noise was acceptably 
small. 

In the RN-DMC process, we can also pick up the sta- 
tistical signal from the walkers that have never crossed 
the nodal surface, and in essence this provides the FN- 
DMC estimator. By monitoring these paths as well, we 
can assess the consistency of the estimators and some- 
what better tune the parameter a for providing better 
RN-DMC signal. 



III. RESULTS 

For benchmark purposes we first calculate perhaps the 
smallest nontrivial system — four atoms. Our result is 
shown in Fig. [T] upper left panel. There is approximately 
10% energy drop when Rcif/rg is reduced from 0.1279 to 
0.003998. We extrapolate i?cfr to zero using a linear fit 
and obtain ^2.2 = 0.212(2). Here and in the rest of the 
paper, the denominator i?froc in the ratio f = E / E^cc is 
evaluated in the same finite volume subject to the same 
boundary conditions as the nominator E. 

The ground-state energy of this small and relatively 
simple system was obtained also by two other numerical 
methods using a lattice formulation of the unitary Fermi 
gas model: the iterative Lanczos diagonalization and the 
auxiliary-field projection Monte Carlo method [55], The 
agreement between our fixed-node DMC results and the 
outcome of these entirely different methods strongly sug- 
gests that our range-extrapolated total energy of the 
four-atom system is very accurate. This conclusion is 
further corroborated by our released-node DMC results 
discussed below. 

Our calculations with 14, 38 and 66 atoms are car- 
ried out analogously to the 4-atom case, and the extrap- 
olated values of f are 0.407(2), 0.409(3) and 0.398(3), 
respectively, as plotted in Fig. [l] In these calculations, 
the smallest effective range is Rcs/fs — 0.003125. It 
should be noted that in the previous DMC calculations, 
the range of Rcs/rs was between 0.1 to 0.2, and within 
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FIG. 2. The pair orbitals and FN-DMC and RN-DMC energies of the 4-atom unitary system with Res/vs = 0.06397. The 
upper row shows the pair orbitals with the lowest (left), intermediate (middle) and optimal (right) accuracy with regard to the 
variational optimization. The lower row shows the corresponding DMC energies as functions of the projection time starting 
from the variational estimate. Note that the resolution of the left and right panels differs by an order of magnitude. The 
vertical dotted lines indicate the instant of the nodal release. 



this range our DMC results agree well with the previously 
obtained values [111 [TH [18] . Reduction of i?eff decreases 
the energy in all cases although the decrease per atom is 
smaller in larger systems. 

To test the quality of the nodal surfaces, we have car- 
ried out released-node calculations for 4-atom and 14- 
atom systems. The RN-DMC calculations for 4 atoms 
were done with Rcff/rg — 0.06397. In a typical released- 
node run the number of walkers was about two million 
so that the error bars were initially very small. In Fig. [2] 
the upper row shows the pair orbital along three distinct 
directions (100, 110 and 111) of the interparticle dis- 
tance vector Ti — Yj . The lower row shows the FN-DMC 
and RN-DMC energies as they evolve with the projection 
time. The plots show convergence of the FN-DMC en- 
ergy followed by the nodal release. This is accomplished 
by switching the guiding wave function from ^'y(R) to 



the bosonic function ^'g'(R) defined in Eq. (11). 



The released-node signal reflects the quality of the 
nodal surface of the trial wave function employed in the 
FN-DMC simulation. We have tested wave functions 
with intentionally varied accuracy by employing subop- 
timal pair orbitals. The plot of the energy evolution in 
the left panel of Fig. [2] shows a clear and pronounced 
drop after the nodal release. As the quality of the pair 
orbital improves, this drop shrinks. For the fully op- 



timized BCS-Jastrow wave function (the right panel in 
Fig. [2]) the energy is reduced by less than 0.002 within 
the longest projection time we have tried. This fact as 
well as comparison with other methods ,33J indicate that 
our BCS wave functions are very accurate in this small 
system and that the fixed-node error is marginal. 

We observe an unexpectedly high sensitivity of the 
nodal quality to the details of the pair orbital at large 
distances. This suggests an explanation for the rela- 
tively slow convergence of the released-node energy: the 
long-range tails of the pair orbital affect the nodal hy- 
persurfaces, although their contribution to the total en- 
ergy is relatively small. One can further deduce that this 
makes the released-node method quite challenging to ap- 
ply since it requires sampling of long distances and the 
corresponding correlations. This is, however, difficult to 
achieve because the diffusive motion of walkers is slow, 
proportional to i^/^, while the growth of the noise is fast, 
proportional to exp(ABi?i) where IS.bf is the difference 
between the bosonic and fermionic ground-state energies. 

The time step At was set to 4 x lO^^r^ in all runs 
and we have verified that the time-step bias of the RN- 
DMC results is negligible. The converged RN-DMC en- 
ergy should not depend on the parameter a in the bosonic 
guiding function "^q provided one would be able to evolve 
the stochastic process with the error bars under control 
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FIG. 3. Evolution of the DMC energies for the 14-atom sys- 
tem with the best optimized BCS-Jastrow wave function. The 
runs are for Rcs/ra = 0.2. No statistically significant energy 
drop is observed after the nodal release that is indicated with 
the vertical dotted line. 



FIG. 4. DMC energies for 14 atoms obtained with the Slater- 
Jastrow wave function. The runs are for Rcs/fs — 0.2. The 
RN-DMC energy drops are significant when compared to the 
RN-DMC signal from the BCS-Jastrow wave function. The 
parameter ^7,7 drops by ~ 0.015 within Ept ~ 0.2 after the 
nodal release. 



until the full convergence. In the same time, the parame- 
ter a crucially affects the growth of the fluctuations with 
the projection time as illustrated in Fig. [2j 

The RN-DMC energy for 14 atoms with Ras/r, = 0.2 
is shown in Fig. [3j The error bars are estimated from 
eight independent runs with two million walkers each. In 
the interval of Ept < 0.2 after the nodal release the RN- 
DMC energy gain appears to be very small and the error 
bars preclude to make any statistically sound estimation 
for longer projection times. The rapid loss of resolution 
is expected since the difference between the bosonic and 
fermionic ground states grows with the number of atoms. 
Again, the RN-DMC signal exhibits little dependence on 
a wc choose. 

In order to make a comparison with a case displaying a 
clear fixed-node bias, we have carried out RN-DMC runs 
using the Slater- Jastrow trial wave function, see Fig. [4] 
Since this wave function has the nodal surface of the non- 
interacting Fermi gas, the nodal surface is strongly dis- 
torted. As a result, we see a very pronounced released- 
node signal. However, within the projection time inter- 
val of Ept < 0.2, the energy drops by only « 0.015 for 
the largest a we tested. This is very small considering 
that the true ground-state energy is at least an order of 
magnitude lower. This again illustrates the challenges 
of efficient application of the released-node method, at 
least for present cases. Comparison of the Slater and 
BCS wave functions shows a significant effect of pairing 
as was demonstrated in earlier studies 

In order to quantify the pairing effects we calculate 
the two-body density matrix which enable us to evaluate 
the condensate fraction. The projected two-body density 
matrix for spin-up and spin-down atoms is defined as 



'I'T(ri,r2) 



(13) 



4F2 /dR$(R)*T(R) 

where N is the total number of atoms and V is the vol- 
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FIG. 5. The two-body density matrix for 66 atoms calculated 
from the FN-DMC mixed estimator. The condensate fraction 
converges to 0.56(1) for Rcs/rs < 0.05. 



ume of the simulation cell. The density matrices have 
been calculated for the fixed-node wave functions and 
hence they correspond to the mixed estimators [25] . Nev- 
ertheless, the mixed-estimator bias is negligible since the 
variational Monte Carlo and DMC estimates of p^^-* co- 
incide within error bars. This is a further evidence of the 
high accuracy of our trial wave functions. 

The condensate fraction can be extracted from the two- 
body density matrix as 



2V^ 



lim p(^)( 



(14) 



The calculated density matrices are shown in Fig. [5] with 
the condensate fraction estimated from the long-range 
limit. The condensate fraction saturates for i?off < 0.5 
at c = 0.56(1). This value is not too far from the results 
obtained previously [TBI I18j . 

To illustrate the character of the nodal surfaces in the 
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FIG. 6. Three-dimensional subsets of the nodal hypersurfaces 
for three types of wave functions and corresponding phases in 
the 14-atom system. The node is obtained by scanning the 
simulation cell with a pair of spin-up and spin-down atoms 
sitting on the top of each other while keeping the rest of the 
atoms at fixed positions (tiny spheres). From the left to the 
right, the columns show the nodal surfaces of the wave func- 
tions corresponding to the free Fermi gas, the unitary limit 
and the BEG side of the crossover. The lower row displays 
the same surfaces rotated by 45 degrees around the z-axis. 



BCS-BEC systems, we present three-dimensional scans 
of the nodes for three wave functions corresponding to 
the following scattering regimes: the free atomic gas with 
no pairing, our best unitary-limit wave function, and the 
wave function with enhanced pairing from the BEC side 
of the BEC-BCS phase diagram {askp = 0.6592). The 
left column of Fig. |6] displays the nodal surface of the free 
atomic Fermi gas. The delocalizcd nature of the system 
is apparent. At the unitary limit, shown in the middle 
column of Fig. |6) the shape of the nodal surface is sig- 
nificantly different as the pairing effects clearly dominate 
and lead to a localized character of the nodes from the 
perspective of a pair of up and down spin atoms. The 
nodes on the BEC side (the right column) do not differ 
much from the unitary limit, except for a slightly more 



pronounced localization. 

IV. CONCLUSIONS 

We have carried out QMC calculations of the zero- 
temperature, spin-unpolarized atomic Fermi gas in the 
unitary limit. We show that the interaction range im- 
pacts the resulting total energies significantly. By ex- 
trapolating the interaction range to zero we obtain the 
ratio E/Efrco for 4, 14, 38 and 66 atoms to be 0.212(2), 
0.407(2), 0.409(3) and 0.398(3), respectively. Our results 
agree well with the previous fixed-node DMC calculations 
when we employ similar simulation parameters such as 
the atom density, the interaction range and the number 
of atoms. From the released-node DMC calculations for 
4 and 14 atoms we have found that the convergence to 
the correct and asymptotically exact ground-state ener- 
gies is unfavorably slow compared to the growth of the 
statistical noise. We were able to identify only small en- 
ergy gains within the simulation times that allowed for 
acceptable signal to noise ratio. We have calculated the 
two-body density matrix and the condensate fraction in 
the limit of zero interaction range, and we have found 
only small changes in these quantities when compared 
with the previous calculations. Our condensate fraction 
from the fixed-node DMC simulations is 0.56(1). 

During the preparation of our manuscript we became 
aware of a similar study where an interaction-range ex- 
trapolation was also performed |34j. The final result for 
the 66-atoni system was ^33,33 = 0.383(2), which is ap- 
proximately 4% lower than ours. We believe that a large 
portion of this difference can be attributed to differences 
in the functional forms of the pairing orbital. Some in- 
fiuence could come also due to the differences between 
the extrapolation methods employed in this work and in 
Ref. [34]. 
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